Convergence Clubs of wage in 34 Indonesian Provinces

#Activate necessary libraries

library(ConvergenceClubs)

package 㤼㸱ConvergenceClubs㤼㸲 was built under R version 4.0.4

library(tidyverse)

Registered S3 methods overwritten by 'dbplyr':
  method         from
  print.tbl_lazy     
  print.tbl_sql      
-- Attaching packages --------------------------------------- tidyverse 1.3.0 --
v ggplot2 3.3.5     v purrr   0.3.4
v tibble  3.1.2     v dplyr   1.0.7
v tidyr   1.1.0     v stringr 1.4.0
v readr   1.3.1     v forcats 0.5.0
package 㤼㸱ggplot2㤼㸲 was built under R version 4.0.5package 㤼㸱tibble㤼㸲 was built under R version 4.0.5package 㤼㸱tidyr㤼㸲 was built under R version 4.0.2package 㤼㸱dplyr㤼㸲 was built under R version 4.0.5-- Conflicts ------------------------------------------ tidyverse_conflicts() --
x dplyr::filter() masks stats::filter()
x dplyr::lag()    masks stats::lag()

library(janitor)

Attaching package: 㤼㸱janitor㤼㸲

The following objects are masked from 㤼㸱package:stats㤼㸲:

    chisq.test, fisher.test

library(tibble)
library(data.table)

package 㤼㸱data.table㤼㸲 was built under R version 4.0.4Registered S3 method overwritten by 'data.table':
  method           from
  print.data.table     
data.table 1.14.0 using 4 threads (see ?getDTthreads).  Latest news: r-datatable.com

Attaching package: 㤼㸱data.table㤼㸲

The following objects are masked from 㤼㸱package:dplyr㤼㸲:

    between, first, last

The following object is masked from 㤼㸱package:purrr㤼㸲:

    transpose

library(ExPanDaR)

Registered S3 methods overwritten by 'htmltools':
  method               from         
  print.html           tools:rstudio
  print.shiny.tag      tools:rstudio
  print.shiny.tag.list tools:rstudio
Registered S3 method overwritten by 'htmlwidgets':
  method           from         
  print.htmlwidget tools:rstudio

#Part I: Nominal Wage using definition “A” of worker #Load the data

library(readxl)
nominal.wage <- read_excel("WS Labour_new.xlsx", 
    sheet = "nominal_wage_A")
nominal.wage

#Select data only from 2008 onwards and remove “Region” column, and remove the observation of Indonesia

nominal.wage <- nominal.wage %>% 
  select(-c(1,3:5)) %>% 
  subset(Province!="National")

#transform the wage series to log

nominal.wage.log <- log(nominal.wage[,-1])
nominal.wage.log

#Remove short-run noise using HP filter

f.nominal.wage.log <- apply(nominal.wage.log, 1, 
function(x){mFilter::hpfilter(x, freq=6.25, type="lambda")$trend} ) 
f.nominal.wage.log <- data.frame(Province = nominal.wage[,1], t(f.nominal.wage.log), stringsAsFactors=FALSE ) 
colnames(f.nominal.wage.log) <- colnames(nominal.wage) 

#Inspect the filtered data

head(f.nominal.wage.log)

#Run the log-t-test

H.nominal <- computeH(f.nominal.wage.log[,-1], quantity = "H")
round(estimateMod(H.nominal, time_trim=0.333, HACmethod = "FQSB"), 3)

    beta  std.err   tvalue   pvalue 
  -1.132    0.011 -100.550    0.000 

#find convergence clubs

clubs.nominal <- findClubs(f.nominal.wage.log, dataCols=2:14, unit_names = 1, refCol=14,
time_trim=0.333, cstar=0, HACmethod = 'FQSB')

summary(clubs.nominal)

Number of convergence clubs: 2
Number of divergent units: 0

        | # of units  | beta     | std.err  | tvalue   | cstar  
-------- ------------- ---------- ---------- ---------- --------
 club1  |   12        |  -0.378  |   0.087  |  -4.316  |   0    
 club2  |   22        |  -0.058  |   0.104  |  -0.558  |   0    

print(clubs.nominal)

================================================= 
club 1
------------------------------------------------- 
Banten, Papua, East Kalimantan, Riau Islands, 
North Kalimantan, Jakarta, West Java, North 
Sulawesi, Central Kalimantan, West Papua, Bali, 
South Sulawesi

beta:     -0.3775 
std.err:   0.0875 
tvalue:   -4.3161 
pvalue:    0 
cstar:     0 

================================================= 
club 2
------------------------------------------------- 
South Kalimantan, Bangka Belitung, North Maluku, 
Maluku, West Sumatera, Riau, Southeast Sulawesi, 
Bengkulu, Aceh, West Kalimantan, Central 
Sulawesi, Gorontalo, North Sumatera, West Nusa 
Tenggara, East Java, Jambi, Lampung, Yogyakarta, 
Central Java, East Nusa Tenggara, South 
Sumatera, West Sulawesi

beta:     -0.0579 
std.err:   0.1037 
tvalue:   -0.5585 
pvalue:    0.2883 
cstar:     0 

plot(clubs.nominal)

plot(clubs.nominal, clubs=NULL, avgTP = TRUE, legend=TRUE)

#Merge clubs

mclubs.nominal <- mergeClubs(clubs.nominal, mergeMethod='PS')
summary(mclubs.nominal)

Number of convergence clubs: 2
Number of divergent units: 0

        | merged clubs  | # of units  | beta     | std.err  | tvalue   
-------- --------------- ------------- ---------- ---------- ----------
 club1  | clubs: 1      |   12        |  -0.378  |   0.087  |  -4.316  
 club2  | clubs: 2      |   22        |  -0.058  |   0.104  |  -0.558  

mclubs.nominal

================================================= 
club 1
------------------------------------------------- 
Banten, Papua, East Kalimantan, Riau Islands, 
North Kalimantan, Jakarta, West Java, North 
Sulawesi, Central Kalimantan, West Papua, Bali, 
South Sulawesi

beta:     -0.3775 
std.err:   0.0875 
tvalue:   -4.3161 
pvalue:    0 

================================================= 
club 2
------------------------------------------------- 
South Kalimantan, Bangka Belitung, North Maluku, 
Maluku, West Sumatera, Riau, Southeast Sulawesi, 
Bengkulu, Aceh, West Kalimantan, Central 
Sulawesi, Gorontalo, North Sumatera, West Nusa 
Tenggara, East Java, Jambi, Lampung, Yogyakarta, 
Central Java, East Nusa Tenggara, South 
Sumatera, West Sulawesi

beta:     -0.0579 
std.err:   0.1037 
tvalue:   -0.5585 
pvalue:    0.2883 

plot(mclubs.nominal)

plot(mclubs.nominal, clubs=NULL, avgTP = TRUE, legend=TRUE)

---

#Part II: Real Wage using definition “A” of worker #Load the data

library(readxl)
real.wage <- read_excel("WS Labour_new.xlsx", 
    sheet = "real_wage_A")
real.wage

#Select data only from 2008 onwards and remove “Region” column, and remove the observation of Indonesia

real.wage <- real.wage %>% 
  select(-c(1,3:5)) %>% 
  subset(Province!="National")

#transform the wage series to log

real.wage.log <- log(real.wage[,-1])
real.wage.log

#Remove short-run noise using HP filter

f.real.wage.log <- apply(real.wage.log, 1, 
function(x){mFilter::hpfilter(x, freq=6.25, type="lambda")$trend} ) 
f.real.wage.log <- data.frame(Province = real.wage[,1], t(f.real.wage.log), stringsAsFactors=FALSE ) 
colnames(f.real.wage.log) <- colnames(real.wage) 

#Inspect the filtered data

head(f.real.wage.log)

#Run the log-t-test

H.real <- computeH(f.real.wage.log[,-1], quantity = "H")
round(estimateMod(H.real, time_trim=0.333, HACmethod = "FQSB"), 3)

   beta std.err  tvalue  pvalue 
 -1.130   0.014 -79.721   0.000 

#find convergence clubs

clubs.real <- findClubs(f.real.wage.log, dataCols=2:14, unit_names = 1, refCol=14,
time_trim=0.333, cstar=0, HACmethod = 'FQSB')

summary(clubs.real)

Number of convergence clubs: 3
Number of divergent units: 0

        | # of units  | beta     | std.err  | tvalue   | cstar  
-------- ------------- ---------- ---------- ---------- --------
 club1  |   3         |   0.113  |   0.232  |   0.486  |   0    
 club2  |   9         |   0.745  |   0.242  |   3.081  |   0    
 club3  |   22        |  -0.014  |   0.108  |  -0.126  |   0    

print(clubs.real)

================================================= 
club 1
------------------------------------------------- 
Jakarta, Riau Islands, Banten

beta:      0.1128 
std.err:   0.232 
tvalue:    0.4863 
pvalue:    0.6866 
cstar:     0 

================================================= 
club 2
------------------------------------------------- 
Papua, East Kalimantan, North Kalimantan, West 
Java, West Papua, North Sulawesi, Bali, Central 
Kalimantan, South Sulawesi

beta:      0.7446 
std.err:   0.2417 
tvalue:    3.0808 
pvalue:    0.999 
cstar:     0 

================================================= 
club 3
------------------------------------------------- 
North Maluku, Riau, Maluku, West Sulawesi, South 
Kalimantan, Bangka Belitung, Southeast Sulawesi, 
West Sumatra, Aceh, Gorontalo, East Java, North 
Sumatra, Bengkulu, Yogyakarta, Central Java, 
South Sumatra, West Nusa Tenggara, Jambi, 
Lampung, Central Sulawesi, West Kalimantan, East 
Nusa Tenggara

beta:     -0.0136 
std.err:   0.1079 
tvalue:   -0.1261 
pvalue:    0.4498 
cstar:     0 

plot(clubs.real)

plot(clubs.real, clubs=NULL, avgTP = TRUE, legend=TRUE)

#Merge clubs

mclubs.real <- mergeClubs(clubs.real, mergeMethod='PS')
summary(mclubs.real)

Number of convergence clubs: 3
Number of divergent units: 0

        | merged clubs  | # of units  | beta     | std.err  | tvalue   
-------- --------------- ------------- ---------- ---------- ----------
 club1  | clubs: 1      |   3         |   0.113  |   0.232  |   0.486  
 club2  | clubs: 2      |   9         |   0.745  |   0.242  |   3.081  
 club3  | clubs: 3      |   22        |  -0.014  |   0.108  |  -0.126  

mclubs.real

================================================= 
club 1
------------------------------------------------- 
Jakarta, Riau Islands, Banten

beta:      0.1128 
std.err:   0.232 
tvalue:    0.4863 
pvalue:    0.6866 

================================================= 
club 2
------------------------------------------------- 
Papua, East Kalimantan, North Kalimantan, West 
Java, West Papua, North Sulawesi, Bali, Central 
Kalimantan, South Sulawesi

beta:      0.7446 
std.err:   0.2417 
tvalue:    3.0808 
pvalue:    0.999 

================================================= 
club 3
------------------------------------------------- 
North Maluku, Riau, Maluku, West Sulawesi, South 
Kalimantan, Bangka Belitung, Southeast Sulawesi, 
West Sumatra, Aceh, Gorontalo, East Java, North 
Sumatra, Bengkulu, Yogyakarta, Central Java, 
South Sumatra, West Nusa Tenggara, Jambi, 
Lampung, Central Sulawesi, West Kalimantan, East 
Nusa Tenggara

beta:     -0.0136 
std.err:   0.1079 
tvalue:   -0.1261 
pvalue:    0.4498 

plot(mclubs.real)

plot(mclubs.real, clubs=NULL, avgTP = TRUE, legend=TRUE)

class(mclubs.nominal)

[1] "convergence.clubs" "list"             

table.mclubs.nom <- map(mclubs.nominal, as.data.table)

Item 3 has 4 rows but longest item has 22; recycled with remainder.

df.mclubs.nom <- rbindlist(table.mclubs.nom, fill = T, idcol = T)
df.mclubs.nom

colnames(df.mclubs.nom)[5] <- "Province"
df.mclubs.nom

nominal.wage.club <- as.data.frame(
  inner_join(f.nominal.wage.log, df.mclubs.nom, by="Province") %>% 
  select(-contains(c("id","model")))
  ) 
nominal.wage.club

#Generate relative value for all clubs

relative <- list()
for(a in 2:14) {
  relative[[a]] <- data.frame(nominal.wage.club$Province, nominal.wage.club[,a]/mean(nominal.wage.club[,a]))
colnames(relative[[a]])[1]<- "Province"
colnames(relative[[a]])[2]<- 
paste("rel",colnames(nominal.wage.club)[a],sep="_")
}

for (x in 2:13) {
  relative[[x+1]]<-left_join(relative[[x]],relative[[x+1]], by="Province")
}

table.relative <- map(relative, as.data.table)
df.relative <- rbindlist(table.relative, fill = T, idcol = T) %>% drop_na()
library(plyr)

package 㤼㸱plyr㤼㸲 was built under R version 4.0.4----------------------------------------------------
You have loaded plyr after dplyr - this is likely to cause problems.
If you need functions from both plyr and dplyr, please load plyr first, then dplyr:
library(plyr); library(dplyr)
----------------------------------------------------

Attaching package: 㤼㸱plyr㤼㸲

The following objects are masked from 㤼㸱package:dplyr㤼㸲:

    arrange, count, desc, failwith, id, mutate,
    rename, summarise, summarize

The following object is masked from 㤼㸱package:purrr㤼㸲:

    compact

df.relative <- join(df.relative,nominal.wage.club,by=c('Province')) %>% 
  select(Province,contains("rel"),clubs)
  
df.relative <- as.data.frame(df.relative)

colnames(df.relative) <- gsub("rel_","",colnames(df.relative))
df.relative

df.relative.long <- df.relative %>%  pivot_longer(-c(Province,clubs), names_to = "Time", values_to="Rel_Wage")
df.relative.long

df.relative.path <- aggregate(Rel_Wage ~ clubs + Time, df.relative.long, mean) %>% 
  arrange(clubs)

path_all <- df.relative.path %>% 
  ggplot(aes(x=Time,y=Rel_Wage, group=clubs, col=clubs)) + geom_line() +
  labs(title = "Transition path of all clubs") +
  theme_bw()
path_all
ggsave("path_all.png", width = 9, height = 5)

#Generate relative value for within club #Create dataframe for each club

club1 <- nominal.wage.club %>% 
  filter(clubs == "club1")

relative.club1 <- list()
for(a in 2:14) {
  relative.club1[[a]] <- data.frame(club1$Province, club1[,a]/mean(club1[,a]))
colnames(relative.club1[[a]])[1]<- "Province"
colnames(relative.club1[[a]])[2]<- 
paste("rel",colnames(club1)[a],sep="_")
}

for (x in 2:13) {
  relative.club1[[x+1]]<-left_join(relative.club1[[x]],relative.club1[[x+1]], by="Province")
}

table.relative.club1 <- map(relative.club1, as.data.table)
df.relative.club1 <- rbindlist(table.relative.club1, fill = T, idcol = T) %>% drop_na()
library(plyr)
df.relative.club1 <- join(df.relative.club1,club1,by=c('Province')) %>% 
  select(Province,contains("rel"),clubs)
  
df.relative.club1

colnames(df.relative.club1) <- gsub("rel_","",colnames(df.relative.club1))
df.relative.club1

df.relative.club1.long <- df.relative.club1 %>%  pivot_longer(-c(Province,clubs), names_to = "Time", values_to="Rel_Wage")
df.relative.club1.long

df.relative.club1.path <- aggregate(Rel_Wage ~ Province + Time, df.relative.club1.long, mean) %>% 
  arrange(Province)

path_club1 <- df.relative.club1.path %>% 
  ggplot(aes(x=Time,y=Rel_Wage, group=Province, col=Province)) + geom_line() +
  labs(title = "Transition path of regions in Club 1") +
  theme_bw() +
   theme(legend.title = element_text(size = 8), 
               legend.text = element_text(size = 8))
path_club1
ggsave("path_club1.png", width = 9, height = 5)

club2 <- nominal.wage.club %>% 
  filter(clubs == "club2")

relative.club2 <- list()
for(a in 2:14) {
  relative.club2[[a]] <- data.frame(club2$Province, club2[,a]/mean(club2[,a]))
colnames(relative.club2[[a]])[1]<- "Province"
colnames(relative.club2[[a]])[2]<- 
paste("rel",colnames(club2)[a],sep="_")
}

for (x in 2:13) {
  relative.club2[[x+1]]<-left_join(relative.club2[[x]],relative.club2[[x+1]], by="Province")
}

table.relative.club2 <- map(relative.club2, as.data.table)
df.relative.club2 <- rbindlist(table.relative.club2, fill = T, idcol = T) %>% drop_na()
library(plyr)
df.relative.club2 <- join(df.relative.club2,club2,by=c('Province')) %>% 
  select(Province,contains("rel"),clubs)
  
df.relative.club2

colnames(df.relative.club2) <- gsub("rel_","",colnames(df.relative.club2))
df.relative.club2

df.relative.club2.long <- df.relative.club2 %>%  pivot_longer(-c(Province,clubs), names_to = "Time", values_to="Rel_Wage")
df.relative.club2.long

df.relative.club2.path <- aggregate(Rel_Wage ~ Province + Time, df.relative.club2.long, mean) %>% 
  arrange(Province)

path_club2 <- df.relative.club2.path %>% 
  ggplot(aes(x=Time,y=Rel_Wage, group=Province, col=Province)) + geom_line() +
  labs(title = "Transition path of regions in Club 2") +
  theme_bw() +
  theme(legend.title = element_text(size = 7), 
               legend.text = element_text(size = 7))
path_club2
ggsave("path_club2.png", width = 9, height = 5)

#Part II: Real Wage using definition “A” of worker #Load the data

library(readxl)
real.wage <- read_excel("WS Labour_new.xlsx", 
    sheet = "real_wage_A")
real.wage

#Regional wage faceted by region

df <- real.wage %>% 
  select(-c(3:5)) %>% 
  subset(Province!="National")
df

#Transform to long format for ploting

df.long <- df %>% 
  pivot_longer(!c(Region, Province), names_to="Year", values_to="Wage")
df.long$Year <- as.factor(df.long$Year)

df.long$Region_f = factor(df.long$Region, levels=c('Sumatra','Java','Bali-Nustra','Kalimantan', 'Sulawesi', 'Maluku-Papua'))

ggplot(df.long, aes(Year, Wage, group=Province)) + geom_line(size=0.3) +
  facet_wrap(~Region_f) +
  theme(
    axis.text.x=element_text(angle=90,hjust=0.5, size=7),
    axis.text.y=element_text(size=7),
    panel.background = element_blank(),
    axis.title=element_text(size=8),
    strip.text = element_text(size=8)
    ) +
  labs(y="Monthly real wage (000 IDR)", size=3)
ggsave("wage.plot.png")

Saving 4.71 x 2.91 in image

#Select data only from 2008 onwards and remove “Region” column, and remove the observation of Indonesia

real.wage <- real.wage %>% 
  select(-c(1,3:5)) %>% 
  subset(Province!="National")

#Transform to long format for ploting

real.wage.long <- real.wage %>% 
  pivot_longer(!Province, names_to="Year", values_to="Wage")

real.wage.long <- as.data.frame(real.wage.long)
real.wage.long$Year <- as.Date(as.factor(real.wage.long$Year), format = "%Y")

real.wage.long$Year <- as.factor(as.numeric(real.wage.long$Year))

prepare_trend_graph(real.wage.long, "Year", c("Wage"))$plot

prepare_quantile_trend_graph(real.wage.long, "Year", c(0.05, 0.25, 0.5, 0.75, 0.95), "Wage")$plot +
  ggplot2::ylab("Real wage (000 IDR)") + scale_x_continuous(breaks = 2008:2020) +
  theme_classic() + guides(color=guide_legend("Quantile")) + theme(axis.text.x=element_text(angle=90, size=10))

#ggsave("wage_quantile.png", width = 9, height = 5)

graph <- prepare_quantile_trend_graph(real.wage.long[c("Year", "Wage")], "Year", c(0.05, 0.25, 0.5, 0.75, 0.95))
graph$plot

#Plot wage of all provinces

graph <- real.wage.long %>%
ggplot(aes(x = Year, y = Wage, group=Province)) +
geom_line(size=0.5, color = "dark gray") +
labs(x = "Time",
y = "Average real wage per month (thousand IDR)", color="") +
theme(text=element_text( family="Palatino")) + theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(),
panel.background = element_blank(), axis.line = element_line(colour = "black")) 
  guides(col = guide_legend(nrow = 8))

$colour
$title
list()
attr(,"class")
[1] "waiver"

$title.position
NULL

$title.theme
NULL

$title.hjust
NULL

$title.vjust
NULL

$label
[1] TRUE

$label.position
NULL

$label.theme
NULL

$label.hjust
NULL

$label.vjust
NULL

$keywidth
NULL

$keyheight
NULL

$direction
NULL

$override.aes
named list()

$nrow
[1] 8

$ncol
NULL

$byrow
[1] FALSE

$reverse
[1] FALSE

$order
[1] 0

$available_aes
[1] "any"

$name
[1] "legend"

attr(,"class")
[1] "guide"  "legend"

attr(,"class")
[1] "guides"

graph
ggsave("wage_all.png", width = 9, height = 5)

#install.packages("ggplotFL", repos="http://flr-project.org/R")
library(ggplotFL)

Loading required package: FLCore
Loading required package: lattice
Loading required package: iterators
FLCore (Version 2.6.15, packaged: 2020-05-26 19:36:47 UTC)

Attaching package: 㤼㸱FLCore㤼㸲

The following objects are masked from 㤼㸱package:dplyr㤼㸲:

    combine, desc, n

The following object is masked from 㤼㸱package:purrr㤼㸲:

    simplify

The following object is masked from 㤼㸱package:tidyr㤼㸲:

    spread

G <- ggplot(real.wage.long, aes(x=Year, y=Wage)) +
  geom_flquantiles(probs=c(0.25, 0.50, 0.75), fill="red", alpha=0.25, show.legend=TRUE) + labs(x = "Time",
y = "Real wage", color="") + theme_classic()
G
ggsave("quantile.png")

Saving 5.26 x 3.25 in image

library(skimr)

package 㤼㸱skimr㤼㸲 was built under R version 4.0.2

desc <- real.wage.long %>% 
  group_by(Province) %>% 
  skim() 

desc <- desc %>% 
  select(!c(4:8))

desc <- desc %>% 
 drop_na() %>% 
  write_excel_csv2("desc_stat.csv")

#Transform the wage series to log

real.wage.log <- log(real.wage[,-1])
real.wage.log

#Remove short-run noise using HP filter

f.real.wage.log <- apply(real.wage.log, 1, 
function(x){mFilter::hpfilter(x, freq=6.25, type="lambda")$trend} ) 
f.real.wage.log <- data.frame(Province = real.wage[,1], t(f.real.wage.log), stringsAsFactors=FALSE ) 
colnames(f.real.wage.log) <- colnames(real.wage) 

#Inspect the filtered data

head(f.real.wage.log)

#Run the log-t-test

H.real <- computeH(f.real.wage.log[,-1], quantity = "H")
round(estimateMod(H.real, time_trim=0.333, HACmethod = "FQSB"), 3)

   beta std.err  tvalue  pvalue 
 -1.130   0.014 -79.721   0.000 

#find convergence clubs

clubs.real <- findClubs(f.real.wage.log, dataCols=2:14, unit_names = 1, refCol=14,
time_trim=0.333, cstar=0, HACmethod = 'FQSB')

summary(clubs.real)

Number of convergence clubs: 3
Number of divergent units: 0

        | # of units  | beta     | std.err  | tvalue   | cstar  
-------- ------------- ---------- ---------- ---------- --------
 club1  |   3         |   0.113  |   0.232  |   0.486  |   0    
 club2  |   9         |   0.745  |   0.242  |   3.081  |   0    
 club3  |   22        |  -0.014  |   0.108  |  -0.126  |   0    

#We find three club convergence fo nominal wage

#Show the membership (provinces) of each club

print(clubs.real)

====================================================== 
club 1
------------------------------------------------------ 
Jakarta, Riau Islands, Banten

beta:      0.1128 
std.err:   0.232 
tvalue:    0.4863 
pvalue:    0.6866 
cstar:     0 

====================================================== 
club 2
------------------------------------------------------ 
Papua, East Kalimantan, North Kalimantan, West Java, 
West Papua, North Sulawesi, Bali, Central Kalimantan, 
South Sulawesi

beta:      0.7446 
std.err:   0.2417 
tvalue:    3.0808 
pvalue:    0.999 
cstar:     0 

====================================================== 
club 3
------------------------------------------------------ 
North Maluku, Riau, Maluku, West Sulawesi, South 
Kalimantan, Bangka Belitung, Southeast Sulawesi, West 
Sumatra, Aceh, Gorontalo, East Java, North Sumatra, 
Bengkulu, Yogyakarta, Central Java, South Sumatra, 
West Nusa Tenggara, Jambi, Lampung, Central Sulawesi, 
West Kalimantan, East Nusa Tenggara

beta:     -0.0136 
std.err:   0.1079 
tvalue:   -0.1261 
pvalue:    0.4498 
cstar:     0 

plot(clubs.real)

plot(clubs.real, clubs=NULL, avgTP = TRUE, legend=TRUE)

#Merge clubs

mclubs.real <- mergeClubs(clubs.real, mergeMethod='PS')
summary(mclubs.real)

Number of convergence clubs: 3
Number of divergent units: 0

        | merged clubs  | # of units  | beta     | std.err  | tvalue   
-------- --------------- ------------- ---------- ---------- ----------
 club1  | clubs: 1      |   3         |   0.113  |   0.232  |   0.486  
 club2  | clubs: 2      |   9         |   0.745  |   0.242  |   3.081  
 club3  | clubs: 3      |   22        |  -0.014  |   0.108  |  -0.126  

mclubs.real

================================================= 
club 1
------------------------------------------------- 
Jakarta, Riau Islands, Banten

beta:      0.1128 
std.err:   0.232 
tvalue:    0.4863 
pvalue:    0.6866 

================================================= 
club 2
------------------------------------------------- 
Papua, East Kalimantan, North Kalimantan, West 
Java, West Papua, North Sulawesi, Bali, Central 
Kalimantan, South Sulawesi

beta:      0.7446 
std.err:   0.2417 
tvalue:    3.0808 
pvalue:    0.999 

================================================= 
club 3
------------------------------------------------- 
North Maluku, Riau, Maluku, West Sulawesi, South 
Kalimantan, Bangka Belitung, Southeast Sulawesi, 
West Sumatra, Aceh, Gorontalo, East Java, North 
Sumatra, Bengkulu, Yogyakarta, Central Java, 
South Sumatra, West Nusa Tenggara, Jambi, 
Lampung, Central Sulawesi, West Kalimantan, East 
Nusa Tenggara

beta:     -0.0136 
std.err:   0.1079 
tvalue:   -0.1261 
pvalue:    0.4498 

plot(mclubs.real)

plot(mclubs.real, clubs=NULL, avgTP = TRUE, legend=TRUE)

#Convert clubs from list object to dataframe

table.mclubs.real <- map(mclubs.real, as.data.table)

Item 2 has 3 rows but longest item has 4; recycled with remainder.Item 4 has 3 rows but longest item has 4; recycled with remainder.Item 3 has 4 rows but longest item has 9; recycled with remainder.Item 3 has 4 rows but longest item has 22; recycled with remainder.

df.mclubs.real <- rbindlist(table.mclubs.real, fill = T, idcol = T)
df.mclubs.real

colnames(df.mclubs.real)[c(1,5)] <- c("Club","Province")
df.mclubs.real

real.wage.club <- as.data.frame(
  inner_join(f.real.wage.log, df.mclubs.real, by="Province") %>% 
  select(-contains(c("clubs","id","model")))
  ) 
real.wage.club

real.wage.club$Club <- gsub("club", "Club ", real.wage.club$Club)
real.wage.club

#Transition path of all clubs

relative.real <- list()
for(a in 2:14) {
  relative.real[[a]] <- data.frame(real.wage.club$Province, real.wage.club[,a]/mean(real.wage.club[,a]))
colnames(relative.real[[a]])[1]<- "Province"
colnames(relative.real[[a]])[2]<- 
paste("rel",colnames(real.wage.club)[a],sep="_")
}

for (x in 2:13) {
  relative.real[[x+1]]<-left_join(relative.real[[x]],relative.real[[x+1]], by="Province")
}

table.relative.real <- map(relative.real, as.data.table)
df.relative.real <- rbindlist(table.relative.real, fill = T, idcol = T) %>% drop_na()
library(plyr)

package 㤼㸱plyr㤼㸲 was built under R version 4.0.4----------------------------------------------------
You have loaded plyr after dplyr - this is likely to cause problems.
If you need functions from both plyr and dplyr, please load plyr first, then dplyr:
library(plyr); library(dplyr)
----------------------------------------------------

Attaching package: 㤼㸱plyr㤼㸲

The following objects are masked from 㤼㸱package:dplyr㤼㸲:

    arrange, count, desc, failwith, id, mutate,
    rename, summarise, summarize

The following object is masked from 㤼㸱package:purrr㤼㸲:

    compact

df.relative.real <- join(df.relative.real,real.wage.club,by=c('Province')) %>% 
  select(Province,contains("rel"),Club)
  
df.relative.real <- as.data.frame(df.relative.real)

colnames(df.relative.real) <- gsub("rel_","",colnames(df.relative.real))
df.relative.real

df.relative.real.long <- df.relative.real %>%  pivot_longer(-c(Province,Club), names_to = "Time", values_to="Rel_Wage")
df.relative.real.long

df.relative.real.path <- aggregate(Rel_Wage ~ Club + Time, df.relative.real.long, mean) %>% 
  arrange(Club)

path_all.real <- df.relative.real.path %>% 
  ggplot(aes(x=Time,y=Rel_Wage, group=Club, col=Club)) + geom_line() +
labs(x = "Time",
y = "Relative transition paths of all clubs", color="") +
theme(text=element_text( family="Palatino")) + theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(),
panel.background = element_blank(), axis.line = element_line(colour = "black")) 
  guides(col = guide_legend(nrow = 8))

$colour
$title
list()
attr(,"class")
[1] "waiver"

$title.position
NULL

$title.theme
NULL

$title.hjust
NULL

$title.vjust
NULL

$label
[1] TRUE

$label.position
NULL

$label.theme
NULL

$label.hjust
NULL

$label.vjust
NULL

$keywidth
NULL

$keyheight
NULL

$direction
NULL

$override.aes
named list()

$nrow
[1] 8

$ncol
NULL

$byrow
[1] FALSE

$reverse
[1] FALSE

$order
[1] 0

$available_aes
[1] "any"

$name
[1] "legend"

attr(,"class")
[1] "guide"  "legend"

attr(,"class")
[1] "guides"

path_all.real
ggsave("path_all.real.png", width = 9, height = 5)

#Transition path of each club #Club 1

club1.real <- real.wage.club %>% 
  filter(Club == "Club 1")

relative.club1.real <- list()
for(a in 2:14) {
  relative.club1.real[[a]] <- data.frame(club1.real$Province, club1.real[,a]/mean(club1.real[,a]))
colnames(relative.club1.real[[a]])[1]<- "Province"
colnames(relative.club1.real[[a]])[2]<- 
paste("rel",colnames(club1.real)[a],sep="_")
}

for (x in 2:13) {
  relative.club1.real[[x+1]]<-left_join(relative.club1.real[[x]],relative.club1.real[[x+1]], by="Province")
}

table.relative.club1.real <- map(relative.club1.real, as.data.table)
df.relative.club1.real <- rbindlist(table.relative.club1.real, fill = T, idcol = T) %>% drop_na()
library(plyr)
df.relative.club1.real <- join(df.relative.club1.real,club1.real,by=c('Province')) %>% 
  select(Province,contains("rel"),Club)
  
df.relative.club1.real

colnames(df.relative.club1.real) <- gsub("rel_","",colnames(df.relative.club1.real))
df.relative.club1.real

df.relative.club1.real.long <- df.relative.club1.real %>%  pivot_longer(-c(Province,Club), names_to = "Time", values_to="Rel_Wage")
df.relative.club1.real.long

df.relative.club1.real.path <- aggregate(Rel_Wage ~ Province + Time, df.relative.club1.real.long, mean) %>% 
  arrange(Province)

a <- path_club1.real <- df.relative.club1.real.path %>% 
  ggplot(aes(x=Time,y=Rel_Wage, group=Province, col=Province)) + geom_line() +
labs(x = "Time",
y = "Relative transition paths of provinces in Club 1", color="") +
theme(text=element_text( family="Palatino")) + theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(),
panel.background = element_blank(), axis.line = element_line(colour = "black")) 
  guides(col = guide_legend(nrow = 8))

$colour
$title
list()
attr(,"class")
[1] "waiver"

$title.position
NULL

$title.theme
NULL

$title.hjust
NULL

$title.vjust
NULL

$label
[1] TRUE

$label.position
NULL

$label.theme
NULL

$label.hjust
NULL

$label.vjust
NULL

$keywidth
NULL

$keyheight
NULL

$direction
NULL

$override.aes
named list()

$nrow
[1] 8

$ncol
NULL

$byrow
[1] FALSE

$reverse
[1] FALSE

$order
[1] 0

$available_aes
[1] "any"

$name
[1] "legend"

attr(,"class")
[1] "guide"  "legend"

attr(,"class")
[1] "guides"

path_club1.real

#ggsave("path_club1.real.png", width = 9, height = 5)

#Club 2

club2.real <- real.wage.club %>% 
  filter(Club == "Club 2")

relative.club2.real <- list()
for(a in 2:14) {
  relative.club2.real[[a]] <- data.frame(club2.real$Province, club2.real[,a]/mean(club2.real[,a]))
colnames(relative.club2.real[[a]])[1]<- "Province"
colnames(relative.club2.real[[a]])[2]<- 
paste("rel",colnames(club2.real)[a],sep="_")
}

for (x in 2:13) {
  relative.club2.real[[x+1]]<-left_join(relative.club2.real[[x]],relative.club2.real[[x+1]], by="Province")
}

table.relative.club2.real <- map(relative.club2.real, as.data.table)
df.relative.club2.real <- rbindlist(table.relative.club2.real, fill = T, idcol = T) %>% drop_na()
library(plyr)
df.relative.club2.real <- join(df.relative.club2.real,club2.real,by=c('Province')) %>% 
  select(Province,contains("rel"),Club)
  
df.relative.club2.real

colnames(df.relative.club2.real) <- gsub("rel_","",colnames(df.relative.club2.real))
df.relative.club2.real

df.relative.club2.real.long <- df.relative.club2.real %>%  pivot_longer(-c(Province,Club), names_to = "Time", values_to="Rel_Wage")
df.relative.club2.real.long

df.relative.club2.real.path <- aggregate(Rel_Wage ~ Province + Time, df.relative.club2.real.long, mean) %>% 
  arrange(Province)

b<-path_club2.real <- df.relative.club2.real.path %>% 
  ggplot(aes(x=Time,y=Rel_Wage, group=Province, col=Province)) + geom_line() +
labs(x = "Time",
y = "Relative transition paths of provinces in Club 2", color="") +
theme(text=element_text( family="Palatino")) + theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(),
panel.background = element_blank(), axis.line = element_line(colour = "black")) 
  guides(col = guide_legend(nrow = 8))

$colour
$title
list()
attr(,"class")
[1] "waiver"

$title.position
NULL

$title.theme
NULL

$title.hjust
NULL

$title.vjust
NULL

$label
[1] TRUE

$label.position
NULL

$label.theme
NULL

$label.hjust
NULL

$label.vjust
NULL

$keywidth
NULL

$keyheight
NULL

$direction
NULL

$override.aes
named list()

$nrow
[1] 8

$ncol
NULL

$byrow
[1] FALSE

$reverse
[1] FALSE

$order
[1] 0

$available_aes
[1] "any"

$name
[1] "legend"

attr(,"class")
[1] "guide"  "legend"

attr(,"class")
[1] "guides"

path_club2.real
ggsave("path_club2.real.png", width = 9, height = 5)

library(plotly)

Attaching package: 㤼㸱plotly㤼㸲

The following objects are masked from 㤼㸱package:plyr㤼㸲:

    arrange, mutate, rename, summarise

The following object is masked from 㤼㸱package:ggplot2㤼㸲:

    last_plot

The following object is masked from 㤼㸱package:stats㤼㸲:

    filter

The following object is masked from 㤼㸱package:graphics㤼㸲:

    layout

ggplotly(b)

`group_by_()` was deprecated in dplyr 0.7.0.
Please use `group_by()` instead.
See vignette('programming') for more help

#Club 3

club3.real <- real.wage.club %>% 
  filter(Club == "Club 3")

relative.club3.real <- list()
for(a in 2:14) {
  relative.club3.real[[a]] <- data.frame(club3.real$Province, club3.real[,a]/mean(club3.real[,a]))
colnames(relative.club3.real[[a]])[1]<- "Province"
colnames(relative.club3.real[[a]])[2]<- 
paste("rel",colnames(club3.real)[a],sep="_")
}

for (x in 2:13) {
  relative.club3.real[[x+1]]<-left_join(relative.club3.real[[x]],relative.club3.real[[x+1]], by="Province")
}

table.relative.club3.real <- map(relative.club3.real, as.data.table)
df.relative.club3.real <- rbindlist(table.relative.club3.real, fill = T, idcol = T) %>% drop_na()
library(plyr)
df.relative.club3.real <- join(df.relative.club3.real,club3.real,by=c('Province')) %>% 
  select(Province,contains("rel"),Club)
  
df.relative.club3.real

colnames(df.relative.club3.real) <- gsub("rel_","",colnames(df.relative.club3.real))
df.relative.club3.real

df.relative.club3.real.long <- df.relative.club3.real %>%  pivot_longer(-c(Province,Club), names_to = "Time", values_to="Rel_Wage")
df.relative.club3.real.long

df.relative.club3.real.path <- aggregate(Rel_Wage ~ Province + Time, df.relative.club3.real.long, mean) %>% 
  arrange(Province)

c <- path_club3.real <- df.relative.club3.real.path %>% 
  ggplot(aes(x=Time,y=Rel_Wage, group=Province, col=Province)) + geom_line() +
labs(x = "Time",
y = "Relative transition paths of provinces in Club 3", color="") +
theme(text=element_text( family="Palatino")) + theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(),
panel.background = element_blank(), axis.line = element_line(colour = "black")) 
  guides(col = guide_legend(nrow = 8))

$colour
$title
list()
attr(,"class")
[1] "waiver"

$title.position
NULL

$title.theme
NULL

$title.hjust
NULL

$title.vjust
NULL

$label
[1] TRUE

$label.position
NULL

$label.theme
NULL

$label.hjust
NULL

$label.vjust
NULL

$keywidth
NULL

$keyheight
NULL

$direction
NULL

$override.aes
named list()

$nrow
[1] 8

$ncol
NULL

$byrow
[1] FALSE

$reverse
[1] FALSE

$order
[1] 0

$available_aes
[1] "any"

$name
[1] "legend"

attr(,"class")
[1] "guide"  "legend"

attr(,"class")
[1] "guides"

path_club3.real

#ggsave("path_club3.real.png", width = 9, height = 5)

ggplotly(c)

Logit model

Load the data

library(readxl)
logit <- read_excel("WS Labour_new.xlsx", 
    sheet = "logit")
logit

Error in gregexpr(calltext, singleline, fixed = TRUE) : 
  regular expression is invalid UTF-8

logit$clubs_real <- factor(logit$clubs_real)
logit$clubs_real08 <- factor(logit$clubs_real08)

logit <- logit %>% 
  mutate(lngrdp = log(grdp))

glimpse(logit)

Rows: 34
Columns: 129
$ Region                 <chr> "SUMATERA", "SUMATERA~
$ Provinsi               <chr> "Aceh", "Sumatera Uta~
$ Province               <chr> "Aceh", "North Sumatr~
$ club                   <dbl> 2, 2, 2, 2, 2, 2, 2, ~
$ clubs                  <dbl> 0, 0, 0, 0, 0, 0, 0, ~
$ club_real              <dbl> 3, 3, 3, 3, 3, 3, 3, ~
$ clubs_real             <fct> 1, 1, 1, 1, 1, 1, 1, ~
$ club_real08            <dbl> 3, 3, 3, 3, 3, 3, 3, ~
$ clubs_real08           <fct> 1, 1, 1, 1, 1, 1, 1, ~
$ manu_gdp               <dbl> 7.362962, 20.058227, ~
$ min_gdp                <dbl> 10.4860421, 1.2372179~
$ export_gdp             <dbl> 6.576492, 29.809383, ~
$ pma_usd_thousand       <dbl> 58750.51, 742857.40, ~
$ pma_share_nat          <dbl> 0.2156826, 2.7687438,~
$ ipm                    <dbl> 71.12000, 72.37333, 7~
$ manu_g                 <dbl> -4.428821, 3.250070, ~
$ min_g                  <dbl> -8.9143761, 7.0445609~
$ export_g               <dbl> 2.241099, 5.021958, 4~
$ pma_usd_g              <dbl> 241.756878, 47.141298~
$ inflation              <dbl> 4.395791, 5.039509, 5~
$ cpi_2008               <dbl> 134.0060, 125.2408, 1~
$ real_wage_2008         <dbl> 1006.5969, 872.7189, ~
$ ump_2008               <dbl> 1000000, 822205, 8000~
$ lbr_prod_2008          <dbl> 64486623, 53444257, 4~
$ yos_2008               <dbl> 8.50, 8.60, 8.26, 8.5~
$ gdp_2008               <dbl> 104597.17, 296095.24,~
$ gdrp                   <dbl> 116485907, 441522073,~
$ gdrp_yoy               <dbl> 1.296245, 5.151837, 5~
$ pmtb_gdrp              <dbl> 31.96295, 29.88666, 2~
$ gdrp_2008              <dbl> 101545237, 331085237,~
$ gdrp_yoy_2008          <dbl> 1.293150, 6.382524, 5~
$ pmtb_gdrp_2008         <dbl> 28.92855, 29.89249, 2~
$ lab_force              <dbl> 2177867.2, 6779890.3,~
$ spa_rwage_knn6         <dbl> 1165.678, 1173.689, 1~
$ spa_rwage_inv          <dbl> 1066.1195, 1075.1873,~
$ tpak                   <dbl> 64.75956, 71.45083, 6~
$ emp_manu               <dbl> 5.058010, 8.072072, 7~
$ emp_agri               <dbl> 44.8070641, 43.037187~
$ emp_cons               <dbl> 6.634741, 5.541747, 5~
$ emp_min                <dbl> 0.7620963, 0.5675358,~
$ emp_trade              <dbl> 16.80033, 20.53636, 2~
$ emp_electricity        <dbl> 0.2737769, 0.2907925,~
$ emp_fin                <dbl> 1.206397, 1.778340, 1~
$ emp_soc_service        <dbl> 20.41898, 14.80034, 1~
$ emp_transcom           <dbl> 4.014083, 5.087083, 4~
$ grdp                   <dbl> 28744.308, 108075.305~
$ grdp_g                 <dbl> 2.975360, 5.576577, 5~
$ geninf                 <dbl> 3.948053, 5.214585, 5~
$ lend_g                 <dbl> 9.778968, 11.768034, ~
$ pcons_gdp              <dbl> 55.37660, 51.73854, 5~
$ inv_gdp                <dbl> 32.20730, 29.78758, 2~
$ totexp_gdp             <dbl> 43.00046, 43.88465, 6~
$ totimp_gdp             <dbl> 52.14954, 35.77006, 6~
$ domexp_gdp             <dbl> 37.990005, 21.226572,~
$ domimp_gdp             <dbl> 50.529657, 23.478050,~
$ forexp_gdp             <dbl> 4.995589, 22.572318, ~
$ forimp_gdp             <dbl> 1.6173898, 12.2851855~
$ agri_gdp               <dbl> 26.73008886, 25.10708~
$ mining_gdp             <dbl> 10.4379411, 1.2370177~
$ manu_gdp_              <dbl> 6.685383, 19.738269, ~
$ constr_gdp             <dbl> 9.113256, 12.237966, ~
$ trade_gdp              <dbl> 16.081582, 19.688706,~
$ finance_gdp            <dbl> 5.778633, 7.987550, 5~
$ gen_serv_gdp           <dbl> 13.982372, 6.691810, ~
$ serv_gdp               <dbl> 19.761005, 14.679360,~
$ dom_open_gdp           <dbl> 0.8851966, 0.4470462,~
$ for_open_gdp           <dbl> 0.06612979, 0.3485750~
$ lend_gdp               <dbl> 27.17989, 38.15352, 3~
$ gdp_labor              <dbl> 57.77792, 71.24029, 6~
$ gdpc_15                <dbl> 23.09289, 31.22043, 2~
$ gdp_sq                 <dbl> 1989.642, 5956.191, 3~
$ pop_sq                 <dbl> 86.15821, 190.77867, ~
$ loan_gr                <dbl> 9.433056, 10.525256, ~
$ spain_knn4_poly        <dbl> 5.014805, 4.714136, 5~
$ spain_cont_poly        <dbl> 5.197999, 4.657841, 5~
$ spain_knn6_cent        <dbl> 5.043715, 4.843269, 4~
$ spain_dist_cent        <dbl> 5.080102, 4.714136, 4~
$ stdev_all_mtm          <dbl> 0.6469880, 0.7405477,~
$ stdev_all_yoy          <dbl> 1.567177, 2.195147, 2~
$ gcons_g                <dbl> 2.301336, 2.860688, 3~
$ gcons_grdp             <dbl> 18.834446, 6.917228, ~
$ grev_g                 <dbl> 9.531354, 16.068736, ~
$ gexp_g                 <dbl> 12.279301, 18.363550,~
$ gdirexp_g              <dbl> 11.66417, 13.77871, 1~
$ gindirexp_g            <dbl> 18.469691, 24.599792,~
$ lbr_prod_g             <dbl> 0.09523408, 3.2211592~
$ trans_cost18           <dbl> 440872.5, 1390255.1, ~
$ std_gdp_g              <dbl> 1.7815710, 0.5159887,~
$ gdp_g_stdev            <dbl> 1.7751960, 0.6697751,~
$ std_pcons_g            <dbl> 0.2628224, 0.4350649,~
$ pcons_g                <dbl> 3.358384, 5.019741, 4~
$ java                   <dbl> 0, 0, 0, 0, 0, 0, 0, ~
$ persistence_1          <dbl> 0.878934, 0.893040, 0~
$ pop_urban              <dbl> 28.11969, 45.21227, 3~
$ poverty                <dbl> 17.931111, 10.430000,~
$ small_micro_emp        <dbl> 7.833679, 4.127258, 7~
$ small_emp              <dbl> 0.8723249, 1.0417838,~
$ micro_emp              <dbl> 7.000823, 3.090453, 6~
$ grdp_gdp_2010          <dbl> 1.4793600, 4.8234093,~
$ grdp_gdp               <dbl> 1.3156809, 4.8910199,~
$ small_output           <dbl> 0.7238740, 2.8516091,~
$ micro_output           <dbl> 1.6259825, 3.0840862,~
$ small_micro_output     <dbl> 2.3498565, 5.9356953,~
$ small_output_gdp       <dbl> 1.3263031, 1.4291156,~
$ micro_output_gdp       <dbl> 4.0170421, 1.9707620,~
$ small_micro_output_gdp <dbl> 5.343345, 3.399878, 7~
$ manu_g_2010_2019       <dbl> -29.40203, 38.02209, ~
$ HHI                    <dbl> 1560.208, 1717.442, 1~
$ unemp                  <dbl> 8.215, 6.366, 6.374, ~
$ wage_hour              <dbl> 13657.8, 12619.6, 137~
$ HHI_avg                <dbl> 0.1578574, 0.1718802,~
$ HHI_g                  <dbl> 0.9210646, 1.0475247,~
$ nom_wage_g             <dbl> 8.357274, 7.442690, 8~
$ nom_wage19_15          <dbl> 35.35608, 32.40019, 3~
$ rel_nom_wage_g         <dbl> 0.79886539, 1.0814938~
$ formal_labor           <dbl> 41.238, 42.062, 36.66~
$ non_formal_job         <dbl> 41.758, 43.062, 49.69~
$ nom_wage_std           <dbl> 1556.2161, 1435.8815,~
$ nom_wage_std_p         <dbl> 1391.9220, 1284.2915,~
$ rel_wage_std           <dbl> 0.8873343, 0.8187211,~
$ unemp_aug_stdev        <dbl> 1.4678877, 0.8854277,~
$ unemp_aug_rel_stdev    <dbl> 1.9673017, 1.1866735,~
$ unemp_std              <dbl> 1.2763652, 0.8086662,~
$ unemp_rel_stdev        <dbl> 1.6937743, 1.0731239,~
$ real_wage_g            <dbl> 7.709113, 6.893801, 7~
$ real_rel_nom_wage_g    <dbl> 0.9565559, 0.8553910,~
$ real_wage_std          <dbl> 1404.6339, 1419.5468,~
$ real_rel_wage_std      <dbl> 0.8940879, 0.9035804,~
$ lngrdp                 <dbl> 10.266195, 11.590584,~

mylogit <- glm(clubs ~ real_wage_2008 + emp_manu + pmtb_gdrp  + tpak + lngrdp, data = logit, family = "binomial")
summary(mylogit)

Call:
glm(formula = clubs ~ real_wage_2008 + emp_manu + pmtb_gdrp + 
    tpak + lngrdp, family = "binomial", data = logit)

Deviance Residuals: 
     Min        1Q    Median        3Q       Max  
-1.71721  -0.12259  -0.01391   0.01037   1.86450  

Coefficients:
                Estimate Std. Error z value Pr(>|z|)
(Intercept)    -62.96071   33.60538  -1.874   0.0610
real_wage_2008   0.03756    0.01823   2.061   0.0393
emp_manu         0.56626    0.34138   1.659   0.0972
pmtb_gdrp        0.72022    0.39324   1.832   0.0670
tpak            -0.23156    0.12841  -1.803   0.0713
lngrdp           1.36255    1.07194   1.271   0.2037
                
(Intercept)    .
real_wage_2008 *
emp_manu       .
pmtb_gdrp      .
tpak           .
lngrdp          
---
Signif. codes:  
0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

(Dispersion parameter for binomial family taken to be 1)

    Null deviance: 44.149  on 33  degrees of freedom
Residual deviance: 13.631  on 28  degrees of freedom
AIC: 25.631

Number of Fisher Scoring iterations: 9